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ABSTRACT 

A three-dimensional hydrodynamical N-body model for the formation of the Galaxy is presented with 
special attention to the formation of the bulge component. Since all previous numerical models for the 
Galaxy formation do not have a proper treatment of the chemical evolution and/or sufficient spatial 
resolutions, we have constructed a detailed model of the chemical and dynamical evolution of the Galaxy 
using our GRAPE-SPH code. Our SPH code includes various physical processes related to the formation 
of stellar systems. Starting with cosmologically motivated initial conditions, we obtain a qualitatively 
similar stellar system to the Galaxy. Then we analyze the chemical and kinematic properties of the bulge 
stars in our model and find qualitative agreement with observational data. The early evolution of our 
model has revealed that most bulge stars form during the sub-galactic merger (merger component of the 
bulge stars). Because of the strong star burst induced by the merger, the metallicity distribution function 
of such stars becomes as wide as observed. We find that another group of the bulge stars forms later in 
the inner region of the disk (non-merger component of the bulge stars). Because of the difference in the 
formation epoch, the main source of iron for this group of stars is different from the merger component. 
Iron in the merger and non-merger components comes mainly from Type II and Type la supernovae, 
respectively. Since a Type la supernova ejects ~ 10 times more iron than a Type II supernova, [Fe/H] 
of the non-merger component tends to be higher than that of the merger component, which widens the 
metallicity distribution function. From these results, we suggest that the Galactic bulge consists of two 
chemically different components; one has formed quickly through the sub-galactic clump merger in the 
proto-galaxy and the other has formed gradually in the inner disk. 

Subject headings: Milky Way: Bulge — galaxy: formation — stars: formation — hydrodynamics 



1. INTRODUCTION 

The formation of bulges, which are the high density cen- 
tral component of spiral galaxies, is a key process in the 
formation of spiral galaxies. There are several scenarios 
for the formation process of the bulges (e.g., Bouwens, 
Cayon, Silk 1999). Such scenarios can be divided into 
the following categories; (1) a primordial collapse (e.g., 
Eggen, Lynden-Bell, & Sandage 1962), (2) merging of 
sub-clumps in a proto-galactic cloud (e.g., Baugh, Cole, 
& Frenk 1996), and (3) a secular evolution of a stellar 
disk (e.g., Norman, Sellwood, & Hasan 1996). Previously, 
the surface brightness of the bulges has been recognized to 
have the i?^/^ profile as in ellipticals, where R denotes the 
distance from the center (e.g., Whitford 1978; Rich 1988). 
However, recent observations have revealed the presence 
of a class of bulges that are well fitted by the exponen- 
tial profile (e.g., Andredakis, Peletier, & Balcells 1995), 
i.e., there exist two classes of bulges according to the sur- 
face brightness. Each formation scenario has advantage 
and disadvantage to explain the diversity of bulges (Wyse 
1999). 

For the formation of the Galactic bulge, Matteucci & 
Brocato (1990) and Matteucci, Romano, & Molaro (1999) 
have constructed one-zone chemical evolution models as 
follows: (1) They assumed that the Galactic bulge was 
formed by a single star-burst. In their best fit model, the 
formation time scale is ~ 0.1 Gyr. (2) One of the predic- 
tions of their calculations is that the element ratio such as 



[0/Fe] remains high even for metal-rich stars ([Fe/H] > 
0) because the formation time scale of the bulge is shorter 
than the lifetime of a Type la supernova progenitor. How- 
ever, the chemical properties predicted by the one-zone 
model depend strongly on their assumption of the dynam- 
ics, i.e., the formation time-scale of the bulge. To constrain 
the formation scenarios, we need to calculate the chemical 
and dynamical evolution together and compare the model 
predictions with both the chemical and kinematic proper- 
ties of the bulge. 

Here, we present our first results of the three- 
dimensional Galactic bulge formation model. Wc model 
the formation of the Galaxy by means of Smoothed Par- 
ticle Hydrodynamics (SPH) method (Lucy 1977; Gingold 
& Monaghan 1977) that incorporates a star formation 
algorithm and a self-consistent treatment of the chemical 
enrichment history of gas. The three-dimensional models 
for the formation of disk galaxies have been investigated 
by many authors (e.g., Katz 1992; Steinmetz & Miiller 
1995; Friedh & Benz 1995; Berczik 1999) and they have 
succeeded in many respects. We follow the method and 
models by the previous authors but use much larger num- 
ber of particles to study the detailed formation and evo- 
lution processes of the Galaxy (Nakasato 2000). 

In Cold Dark Matter (CDM) cosmology, the early stage 
of the formation of galaxies involves the progressive merg- 
ing of sub-galactic clumps (e.g., Baugh, Cole, & Frenk 
1996). To properly model the dynamics of the merger 
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history, wc adopt the three-dimensional SPH method. 
The chemodynamical approach by Samland, Hensler, & 
Theis (1997), which is a mesh based and two-dimensional 
method, is not suitable to model the early merging history 
of the proto-galaxy. 

Our SPH code uses the GRAPE (Sugimoto et al. 1990) 
and the Remote-GRAPE system (Nakasato, Mori, & 
Nomoto 1997) to compute the gravitational interaction 
between particles (GRAPE-SPH code). The adoption of 
the Rcmotc-GRAPE system doubles the performance of 
our code since we can make the SPH calculation and N- 
body summation in parallel (Nakasato, Mori, & Nomoto 
1997). This enhancement of performance enables us to 
use seven times larger number of particles than Katz 
(1992), Steinmetz & Miiller (1995), and Raiteri, Villata, & 
Navarro (1996) so that the spatial resolution of our model 
is almost two times higher than their models. As a result, 
we can investigate the early merging history in much more 
detail during the Galaxy formation. We note that other 
authors (e.g., Navarro & Steinmetz 2000; Sommer-Larsen 
& Dolgov 2001) have published the comparable or even 
higher resolution SPH galaxy formation models than the 
present work; however, these authors have mainly con- 
centrated on the kinematic properties of the models. In 
contrast, in the present work, we have made attention to 
both the chemical and kinematic properties. 

Our study is the first attempt to model the formation 
and chemical evolution of the Galactic bulge using the 
self-consistent three-dimensional SPH method. Using the 
chemical and dynamical SPH code, we simulate the evo- 
lution of the Galaxy starting from a plausible cosmologi- 
cal initial condition. The obtained model well reproduces 
the important chemical and kinematic properties of the 
Galaxy such as the formation of the bulge, disk and halo 
as seen in Figure 5. By analyzing the result of this model, 
we construct a plausible scenario for the Galactic bulge 
formation. 

The plan for the present paper is as follows: In section 2, 
we briefly describe our SPH method. Section 3 describes 
the procedure to construct initial conditions. We present 
a summary of results obtained with our model in section 4. 
By analyzing the result of our Galaxy model, we describe 
the formation history of the bulge component in section 5, 
and compare the chemical and kinematic properties of the 
bulge with recent observational results. Finally, section 6 
is devoted to conclusions. 

2. THE METHOD 

We adopt the SPH method to construct self-consistent 

three-dimensional dynamical and chemical models for the 
formation of stellar systems. The SPH method has been 
applied to many astrophysical problems. Because of its 
Lagrangian nature, it is suitable to systems that have 
large density contrasts, e.g., the formation of galaxies (e.g., 
Evrard 1988; Hernquist & Katz 1989; Katz 1992; Stein- 
metz & Miiller 1994), the evolution of galaxies (e.g., 
Friedli & Benz 1995; Patsis & Athanassoula 2000), the 
cosmological simulations (e.g., Navarro & White 1994; 
Yoshikawa, Jing, & Suto 2000), and a cloud-cloud collision 
(e.g., Lattanzio et al. 1985; Habe & Ohta 1992). Various 
codes have been developed to combine SPH and N-body 
systems. In these codes, gravitational forces are calculated 



with various methods such as direct summations, Particle- 
Particle/Particle-Mesh methods (e.g., Evrard 1988), Tree 
methods (e.g., Hernquist & Katz 1989; Benz et al. 1990), 
and the method to use the special purpose computer 
GRAPE (e.g., Umemura et al. 1993; Steinmetz 1996). 

To simulate the formation of a stellar system from gas, 
we use our GRAPE-SPH code using the Remote-GRAPE 
library (Nakasato, Mori, & Nomoto 1997). The SPH for- 
mulation that we use is the same as Navarro & White 
(1993). Wc use a spatially variable smoothing length and 
integrate equations of motion with a second order Runge- 
Kutta method as described in Navarro & White (1993). 
To simulate the formation of stellar systems from gas, our 
GRAPE-SPH code includes various physical processes re- 
lated to the formation of stellar systems, e.g., cooling, star 
formation, and feedback by stars as follows. 

2.1. Cooing and star formation 

We adopt the metallicity dependent cooling functions, 
which are computed by MAPPINGS HI package by R.S. 
Sutherland (MAPPINGS HI is the successor of MAP- 
PINGS II described in Sutherland & Dopita 1993), to 
solve the energy equation of gas particles (see Figure 1 of 
Nakasato, Mori & Nomoto (2000) for the actually used 
cooling functions) 

For the star formation algorithm, we adopt the star for- 
mation recipe as used in the usual SPH codes (Katz 1992). 
Hereafter, "STAR" means "star particle", which is not an 
individual star but an association of many stars. In this 
recipe, we convert a fraction of a gas particle into a STAR 
when four conditions are satisfied. The first three condi- 
tions are the physical conditions expressed as 

(V-t;), <0, (1) 

<cool < td, (2) 
td *i isound) (3) 

where fcooh ) and tgound are the cooling time, dynami- 
cal time, and sound crossing time of the gas particle, re- 
spectively (sec Nakasato, Mori & Nomoto (2000) for the 
detailed time scale definition). 

During each star formation interval {St = 2 Myr in the 
present work), we check whether the three conditions are 
satisfied for each gas particle. If the gas particle satisfies 
these three conditions, we then check the forth probability 
condition. Assuming star formation process is a Poisson 
process (Katz 1992), we can write the probability for star 
formation within 5t as 

P = 1.0 - exp(-(5t/tstarform), (4) 

where fstarform is the local star formation time scale. In 
our SPH code, we define istarform = t^/C, where C is a 
parameter to control our star formation recipe. If P is 
larger than a random number (0 — 1), we create a new 
STAR from a fraction of the gas particle (star formation 
occurs!). Introducing this probability condition, we can ef- 
fectively control the threshold density for star formation 
by changing the parameter C. We compute P as a func- 
tion of the density (p) for different values of C as shown 
in Figure 1. Here, we assume St = 2 Myr. Note that 
the timestep used in the time integration is always smaller 
than St. From Figure 1, we can see effects of C as follows; 
larger (smaller) C means lower (higher) threshold density. 
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To demonstrate the dependence on the star formation 
parameter C, we evolve a spherical region of 10^^ Mq, 
which consists of 10 % gas and 90 % dark matter, for var- 
ious C. Initially, we set up the sphere in rigid rotation 

and outward expansion. And the angular speed of rota- 
tion corresponds to the spin parameter A ~ 0.05, where 

A = Qj^^5/2 1 with L being the angular momentum, E the 
total energy, M the mass, and G the gravitational constant 
(Padmanabhan 1993). This initial condition is the same 
as used in Katz (1992) and Steinmetz & Miiller (1995). 
We evolve the sphere for three different star formation pa- 
rameters C = 0.1, 0.5 and 1.0. In Figure 2, we present 
the star formation history for C = 0.1 and C = 1.0. For 
C = 1.0, the star formation rates (SFR) increase more 
rapidly than C = 0.1 since the threshold density for star 
formation is lower. These STARs which formed early (i.e., 
very old STARs) in the C = 1.0 model have large verti- 
cal velocities. If we compare the projected edge-on surface 
density of STARs, the C = 1.0 model clearly shows more 
extended distribution then the C = 0.1 model. In Table 
2, some kinematic properties at t = 5 Gyr are compared 
for different C. There is a clear correlation between C and 
the ratio between the X and Z axis velocity-dispersion (or 
the spin parameter). From these results, we conclude that 
the C = 0.1 and C = 1.0 galaxies are similar to late type 
and early type galaxies, respectively. To summarize, if we 
use this probability condition, various values of C lead to 
various global star formation histories. Hereafter, we set 
C = 0.1 throughout the present paper. 

The mass of a newly formed STAR (Mgtar) is calculated 

as 

1 

1 



star 



1 + 0.5 



St 



^starfor] 



(5) 



where h and p are the smoothing length and the den- 
sity of the gas particle, respectively (Nakasato, Mori & 
Nomoto 2000). In the present calculation, the typical 
mass of STARs is ~ 10^ - 10^ Mq. 

2.2. Feedback from stars 

Formed stars eject energy and mass as stellar winds and 
supernova explosions (feedback from stars). As a result, 
stars heat up, accelerate, and enrich the circumstellar and 
interstellar medium. With current computing resources, 
it is not yet feasible to accurately include the energy, mo- 
mentum, and mass feedback from stars in our SPH code. 
This is because the length resolution of the present cal- 
culation (~ 0.5 kpc) is much longer than a typical size 
of supernova remnants (< 100 pc). Thus, in the present 
paper, we adopt a phcnomcnological method as adopted 
by several authors (e.g., Katz 1992) to mimic real feed- 
back processes; i.e., we distribute the energy and mass 
ejected by a STAR to the neighbor gas particles. In our 
SPH code, we incorporate the effect of stellar winds and 
Type II and Type la supernova explosions as the energy 
and mass feedback processes. We compute the mass and 
energy ejection rates based on a simple model explained 
below. In this paper, we distribute the feedback energy as 
purely thermal energy. This treatment is a simple and a 
zero-th order approximation to the real nature and other 
authors have tried to refine the treatment of feedback (e.g., 
Thacker & Couchman 2000; Springel & Hernquist 2002). 



Also in this paper, wc follow the chemical evolution of total 
metal (Z), iron (Fe) and oxygen (O). The chemical evolu- 
tion is coupled with the thermal and dynamical evolution, 
i.e., the ejected metal affects the evolution of gas through 
the mctallicity dependent cooling rates. To summarize, 
our GRAPE-SPH code is a most up-to-date SPH imple- 
mentation for the formation of stellar systems and used 
to successfully model the formation of globular clusters 
(Nakasato, Mori & Nomoto 2000). 

2.3. Energy ejection 
The energy ejection rate per STAR is given as 

-Eojcct = esW-RsW + esNII-RsNII + esNIa-RsNIa, (6) 

where esw is the total ejected energy by stellar winds dur- 
ing the stellar life time and esNii and esNia are the energies 
ejected by a Type II and la supernova explosion, respec- 
tively. The rate i?sw is the number of stars that expel 
their envelopes at the current epoch per unit time and 
i?SNii and iisNia are the rates of Type II and Type la 
supernovae, respectively. We define Rsw and Rsnu as 
follows 

6{m)dm 



Rsw 



-RSNII 



T(Mm,) 



/■M„a 

/ (j){m)dm 
r(M^s) - T(M„a) ' 



(7) 



(8) 



where 4>{m) is the initial mass function (IMF), namely 
(j){m)dm gives the number of stars in the mass range of 
(m, m + dm), and t(to) is the stellar lifetime as a bmction 
of stellar mass (David, Forman & Jones 1990). In the 
present study, we assume the power-law type IMF as 

0(m) = Am-2•3^ (9) 

where A is a constant. For the upper and lower mass lim- 
its in the IMF, M^p = 120 Mq and Mio = 0.05 Mq are 
assumed. In Eqs. (7) and (8), M^a (= 50 Mq) and M^s 
(= 8 Mq) are the upper and lower mass limits of the stars 
that explode as Type II supernovae. 

For the Type la supernova rate, wc adopt the formu- 
lation given in the galactic chemical evolution model by 
Kobayashi et al. (1998, 2000). They adopt the single- 
degenerate scenario (Nomoto et al. 1994) for the Type la 
supernova progenitor, which is that a white dwarf (WD) 
in a close binary undergoes a thermonuclear explosion 
when the companion star evolves off the main-sequence 
and transfers a large enough amount of mass over to the 
WD. In their model, the Type la supernova rate at an 
epoch t is given as 



■RsNIa(i) = C'sNIa 



M{t+At) 



(j}{m)dm 



M(L) 



At 



(10) 



where M(t) is the mass of the companion star whose main- 
sequence lifetime is t and CsNia is a constant to be cali- 
brated by the observational constraints. From the evo- 
lution model of the progenitor stars, Hachisu, Kato & 
Nomoto (1996, 1999) obtained the mass ranges of the 
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companion stars of the WDs that become Type la super- 
novae as 

Msec = [1.8,2.6], [0.9,1.5] Me. (11) 

We note that these mass ranges weakly depend on the 
metalhcity of the star. In the present paper, we don't 
adopt the mctaUicity dependence. Figure 3 shows i?sNia 
as a function of time, where we set CsNia = l-O. In the 
present paper, we set CsNia = 4.0 x lO""* from the result 
of test calculations. Note that other authors (Raiteri, Vil- 
lata, & Navarro 1996; Carraro, Lia & Chiosi 1998) adopt 
different progenitor models by Greggio & Renzini (1983) 
and Matteucci & Greggio (1986). 

For the supernova energy, we assume that csnii = 



eSNIa 



10 erg. For the stellar wind energy, esw is 
estimated to be 0.2x10^^ erg for solar metallicity stars 
from the observational data of OB associations (Abbot 
1982). The chemical abundance of a massive star signif- 
icantly affects esw (Leitherer, Robert, & Drissen 1992), 
so that we use metallicity dependent esw a-s esw = 
0.2esNii(.^/.^o)°'*, where Z is the mass fraction of heavy 
elements in the STAR. 

2.4. Mass ejection 

In our code, the mass ejection in stellar winds from mas- 
sive stars (m > Mms) is combined with the mass ejection 
due to Type II supernovae. The mass ejection in stellar 

winds from low mass stars (m < Mms) is treated sepa- 
rately. Thus, the mass ejection rate per STAR is written 
as 

-Reject = ?TISNII-RSNII + ?TlsWm-RsWm + WSNIa-RsNIa, (12) 

where msNii is the average ejecta mass of Type II super- 
novae, and mswm is the average mass ejected in stellar 
winds from the low mass stars (corresponding to plane- 
tary nebulae), and msNia is the ejecta mass of a Type la 
supernova. i?swm is the number of stars per unit time 
expelling their envelopes at a certain epoch. msNii and 
wswm are defined as 

m0(m)dm 

mNS, (13) 



TOSNII 



(j){m)dm 



mswm 



Jmu 



m^{m)dm 



J Mm, 



(14) 



{m)dm 



Here ttins is the mass locked up in a neutron star and 
mwD is the mass that is locked up in a WD. M\i (= 1 
Mq) is the mass of the star whose life time nearly equals 
to the Hubble time. We set mNS = 1-4 Mq and towd = 
1.4 Mq and assume that a Type la supernova is produced 
by a Chandrasckhar mass WD so that we set mgNia =1-4 
Mq. We define i?swm as 



-RsWm = 



Mil 



(15) 



t(Mu) - t{M^^) 

The fraction of heavy elements in msNii and msNia is com- 
puted using the nucleosynthesis yield of Type II and la 



supernovae (Tsujimoto et al. 1996; Nomoto et al. 1997). 
We compute the chemical evolution of total metal (Z), iron 
(Fe) and oxygen (O) in our code. The metallicity yield we 
used is tabulated in Table 1. 

2.5. Summary on Feedback 

The feedback phase is divided into three phases; a stel- 
lar wind phase, a Type II supernova phase, and a Type la 
supernova phase. 

• The stellar wind phase continues for r(Afnia), dur- 
ing which only the energy ejection from STARs is 
included; the ejected mass is included in the Type 
II supernova phase for simplicity. 

• The Type II supernova phase begins at t = r(Afina) 
and ends aX t ~ T(Mms)- During this phase, the 
mass ejection rate is the sum of the contributions 
by stellar winds of massive stars and Type II super- 
novae. 

• The Type la supernova phase begins &it = t{M^s)- 
During this phase, both the energy ejection and 
mass ejection rates are the sum of the contribu- 
tions by Type la supernovae and stellar winds of 
low mass stars. 

We present the schematic view of the feedback processes 
in Figure 4. 

With the present algorithm and the adopted parame- 
ters, a STAR of IQS Mq produces - 5.5 x 10'"^ Type II 
supernova explosions and 1.4 x 10^ Type la supernova 
explosions during 13 Gyr of the evolution. Also ^ 2.2 x lO'' 
Mq gas is ejected during the evolution. The ejected gas 
contains ~ 1.6 x 10® Mq of heavy elements, which include 
~ 1.5 X 10''' Mq of iron and 1.0 x 10® Mq of oxygen. 

As a final note, the thermal energy, gas, and heavy el- 
ements from stellar winds and supernovae are smoothed 
over neighbor particles of the STAR within a neighbor ra- 
dius oi Rf (a feedback radius). i?j is a parameter set to 
0.5 kpc, which is equal to the adopted gravitational soften- 
ing length of STAR particles. Results of a test calculation 
with Rf = 0.8 kpc are very similar to those of the present 
paper. 

3. INITIAL MODEL 

Following the previous work (e.g., Katz 1992). we inves- 
tigate the evolution of the spherical 3 a top hat over-dense 
region of a mixture of dark matter and gas with our SPH 
code. We neglect the matter outside the spherical regions, 
thus neglecting the late in-fall of the external matter and 
the external gravitational field. Throughout the paper, we 
set the Hubble constant Hq = 50 km s~^ Mpc~^, 0=1, 
and A = 0. We set the initial redshift (z) of the model 
to be ^ 25 and the initial radius of the sphere to be 55 
kpc. Thus, the comoving radius of the spherical region is 
~ 1.4 Mpc and the contained mass is ~ 10^^ Mq. We set 
the region in rigid rotation to provide a sufRcient angular 
momentum with the spin parameter A ^ 0.1. We also add 
the corresponding Hubble velocity to the velocity field of 
the sphere since we integrate the model with the physical 
coordinate. 
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To generate a 3 cr over-dense region, we use the path 
integral method (Bertschinger 1987). We note that there 
are further two methods to set up the over-dense region: 
(1) To seek the desired over-dense region by sampling dif- 
ferent random field realizations, which is used by the cos- 
mological hydrodynamical simulations of formation of X- 
ray clusters (Anninos & Norman 1996). (2) To select the 
desired halo from the results of the large scale cosmological 
N-body simulations and use it as an initial condition for 
a cosmological galaxy formation model (Navarro & White 
1994). The path integral method is the fastest method 
among other methods. And it is the easiest way to gen- 
erate a several dozen desired over-dense regions, i.e. the 
over-dense region of 10^^ Mq in the present study. The 
weakness of this method is that we have to use the isolated 
boundary condition. Namely, we neglect the external tidal 
field and artificially add the angular momentum to the ini- 
tial velocity field. On the other hand, with the method (2), 
it is natural to use the surrounding particles as the origin 
of external tidal- field. However, with the method (2), it 
is difficult to find several over-dense regions with similar 
properties since each region is very different like the real 
universe. We think our way to generate initial conditions 
is an efficient way for our purpose to construct the model 
for our Galaxy. 

Specifically, initial models are constructed as follows: 
(1) We generate a spherical top hat 3 a over-dense region 
by the path integral method (Bertschinger 1987). We use 
the COSMICS package (Bertschinger 1995) and gener- 
ate 50 realizations. (2) We follow the non-linear evolution 
of the dark matter particles in each spherical region with 
the same code. Note that we only use the N-body part 
of our GRAPE-SPH code in these simulations. (3) At a 
certain epoch {z ~ 2.5), we examine the properties of dark 
matter halos and select an appropriate realization for the 
hydrodynamical model. We select the model with a single 
halo since if two dominant halos exist, the gas disk will 
be destroyed by a major merger event (Barnes & Hern- 
quist 1992). Among 50 realizations, 39 realizations are 
single halo models and others are multiple halo models. 
From 39 realizations, we select candidates for the hydro- 
dynamical model by comparing the model density profile 
with the density profile of the Galactic halo. (4) Once we 
select several models, we restart the hydrodynamical N- 
body simulation from ^ ~ 25 with our GRAPE-SPH code. 
After we have done full hydrodynamical calculations for 
three single-halo models (and two multiple-halo models), 
we select the best-fit model, which is presented in this pa- 
per, by comparing the disk scale length with the Galactic 
scale length. All results shown in the following sections 
are obtained with this best-fit model. 

We have computed a similar kind of models with several 
different configurations and found that the performance 
mainly depends on the speed of the GRAPE system: (a) 
For our old configuration using four GRAPE-3A boards, 
a similar calculation to the present work takes more than 
200 hours up to ^; ~ 0.9 (f = 5 Gyr). (b) For our recent 
configuration using one GRAPE-5 board, a similar calcu- 
lation takes only ^ 150 hours up to z = (f = 13 Gyr). 
Using this recent configuration, a full hydrodynamical N- 
body computation for several dozen realizations are under 
way and results will be presented elsewhere. 



4. FORMATION OF THE GALAXY 

The initial number of gas particles and dark matter par- 
ticles are ^ 27, 000 and ^ 27, 000, respectively. Thus, with 
the total mass of 10^^ Mq, the initial masses of the gas 
particles and dark matter particles are 4 x 10^ Mq and 
3.6 X lO'' Mq, respectively. We set the gravitational soft- 
ening length for the gas and dark matter particles to be 
0.5 kpc and 1.0 kpc, respectively. The overall evolution is 
very similar to the previously reported work (Steinmetz & 
Miiller 1995; Berczik 1999). During the expansion due to 
the Hubble flow, small scale structures grow and dense re- 
gions collapse. In such high density regions (clumps), star 
formation occurs as a result of the efficient radiative cool- 
ing. These small clumps gradually merge to produce larger 
clumps. At f ~ 0.3 - 0.5 Gyr, these clumps merge and the 
resulted clump becomes the primary halo that eventually 
becomes a spiral galaxy. Around this time (z ~ 7), the star 
formation rates (SFR) reaches a maximum value. At f ~ 
1 Gyr (2 ~ 4), the dark matter is almost virialized and the 
evolution of dark matter particles becomes a quasi-static 
state while the gas partic;les are being settled into a thin 
disk. Finally, the prominent gas disk forms in the center 
of the primary halo at t ~ 3 Gyr {z 1.6). Afterward, the 
most star formation occurs in this gas disk and the overall 
evolution becomes a quasi-static state. 

After 5 Gyr of evolution (up to z ~ 0.9), over 80,000 
STARs have formed and the number of gas particles have 
decreased to ~ 13, 000. At this stage, we categorize STARs 
according to the chemical properties and the formation 
epoch (tform) of STARs following the observational charac- 
teristics of the three components of the Galaxy as defined 
in Mihalas & Binney (1981). 

Specifically, we categorize STARs by the following con- 
ditions and show the projected stellar densities in Figure 
5: 

• left panel : Metal-poor (log(Z/Z0) < -3) STARs, 

• center panel : Old (fform < 1 Gyr) and metal-rich 
(log(Z/Z0) > 0) STARs, 

• right panel : Young (tform > 4 Gyr) STARs. 

Here, Z is the mass fraction of heavy elements in STARs. 
Metal-poor STARs show an extended distribution since 
these STARs formed at an early epoch and at high lati- 
tude (the left panel). The old and metal-rich STARs are 
concentrated in the galactic center (the center panel) . The 
qualitative properties of the system are not so changed af- 
ter ~ 1 as reported by Steinmetz & Miiller (1995). 
Since we neglect the late in-fall of surrounding material, 
the most recent star formation mainly occurs in the disk. 
As a result, the young STARs are located near the galactic 
plane (the right panel). These results show that the calcu- 
lated stellar system looks very similar to the Galaxy. To 
be more qualitatively, we plot the face-on surface density 
profile of STARs in Figure 7. The solid crosses show our re- 
sult. We fit the model with the following function; p{r) = 
/9diskexp(-r/i?disk) + /0buigeexp(-7.67((r/J?buige)°'^^ - 1)), 
where the scale radiuses are i?disk = 4.21 kpc and i?buigo = 
0.876 kpc. In the next section, we analyze the properties 
of STARs near the galactic center. 

At ^ = 0, the total number of particles has increased 
to 150,000 where the numbers of gas, dark matter, and 
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STAR particles arc - 11,000, 27,000, and 112,000, respec- 
tively. When wc compare the projected density profile at 
z ~ 0.9 (Figure 5) and at ^ = 0, the 3 components struc- 
ture looks very similar except the thickness of the disk 
STARs. Most star formation occurs quasi-statically near 
the disk plane and the SFR is ~ 1 - 2 Mq yr'^ 

5. FORMATION HISTORY OF THE BULGE COMPONENT 

The one-zone chemical evolution model (e.g., Matteucci, 
Romano, & Molaro 1999) predicted the approximate 
chemical properties of the Galactic bulge. To construct 
a self-consistent model of the bulge formation, however, 
a chemical and dynamical model is required. In this sec- 
tion, we use our high resolution model for the formation 
and evolution of the Galaxy to study the formation of the 
Galactic bulge. 

5.1. Dynamical history and kinematics of the bulge stars 

Wc select the STARs located near the galactic center 
{R < 2 kpc) and define these STARs as the bulge stars. 
We define the center of the gravity of all STARs as the 
galactic center. To find the galactic center, we need several 
times of a iteration to compute the center of the gravity. 
In the i-th iteration, we compute the center of the gravity 
of STARs that is located within Ri kpc from the origin of 
the coordinate and convert the origin with that point. We 
set Ri = 20 kpc and R2, Rs, ■■■ = 5 kpc. Typically, we only 
need three or four iterations. For example, at z ^ 0.9, the 
mass and number of the bulge stars are 1.8 x lO-'^" Mq and 
16,800, respectively. Ever after this time, we can select the 
bulge stars with the same definition. In this paper, how- 
ever, we compare our model results at this time with some 
observational data. The reason is that in SPH models in- 
cluding star formation recipes, the resolution of the model 
becomes worse and worse during the evolution since the 
number of gas particles decreases. With our chemical en- 
richment model, we need to have a large enough number 
of gas particles around a STAR particle to properly model 
the chemical evolution of the interstellar matter; further, 
the gas particles show a more extended distribution than 
the star particles and this distribution makes the effective 
resolution of the chemical enrichment description worse. 
At 2 ~ 0.9, the number of gas particles becomes 13,300, 
which is a half of the initial number of gas particles and 
sufficiently large. 

Figure 6 shows the SFR as a function of time for the 
bulge stars (solid line) and all stars (dotted line) up to 
z ~ 0.9. In this figure, we normalize the SFR by the total 
mass of each component at z ~ 0.9. More than 60 % of 
the bulge stars up to this time formed during the first 0.5 
Gyr of the evolution. With our model, we can obtain the 
detailed formation history of the bulge stars. In Figure 8, 
we present snap shots for projected STAR positions in the 
region where star formation is most violent for the first 
0.5 Gyr. We note that the merging of two sub-galactic 
clumps occurs during t = 0.4 - 0.5 Gyr. This merging of 
the clumps causes the star burst of the short time scale 
shown in Figure 6. Such a phenomenon is different from 
the simple star burst assumed in the one-zone chemical 
evolution model (Matteucci & Brocato 1990). 

With our chemical and dynamical model, we can also ob- 
tain the kinematics of the bulge stars. In Figure 9, we plot 



the average velocity (upper line) and the velocity disper- 
sion (lower line) of these bulge stars. The bulge stars are 
substantially rotationally supported with the rotational 
speeds of ~ 160 - 230 km s~^. These speeds are smaller 
than the rotational speeds of the disk stars (^ 270 km s^^ 
in the present model), which means that the angular mo- 
mentum of the bulge is lost during the merging process. 
Figure 10 shows the x-direction velocity dispersion for the 
bulge stars as a function of [Fe/H]. Here, our numerical re- 
sults are shown by the solid triangles and solid squares and 
compared with observational data (Minniti 1996) shown 
by the doted circles (note that the observational data is the 
line of sight velocity dispersion). The solid squares are the 
STARs with tform < 1 Gyr (old STARs) and the solid tri- 
angles are the STARs with t[orm < 5 Gyr (young STARs). 
The old and young STARs show a different trend, and the 
result for old STARs is more consistent with the observa- 
tional data. Our model (especially old STARs) qualita- 
tively reproduces the kinematic properties of the Galactic 
bulge. 

5.2. Chemical properties of bulge stars 

Important information to constrain the formation of the 
Galactic bulge concerns the chemical properties of the 
bulge stars. In this subsection, we analyze the chemical 
properties of the bulge stars. 

The observed mctallicity distribution function of K- 
giant stars (McWilliam & Rich 1994) is characterized by 
a wide (—1.2 < [Fe/H] < 0.9) distribution as shown by 
the dashed line in Figure 11. The distribution function of 
the bulge stars in our model (the solid line in Figure 11) 
is as wide as —2.0 < [Fe/H] < 1.0, being in good agree- 
ment with the observational data. As noted in the previ- 
ous subsection, about 60 % of the bulge stars form in the 
strong star-burst caused by merging of sub-clumps (here- 
after these stars are called a merger component). This 
merger-driven star-burst is partly responsible for this wide 
distribution. On the other hand, other 40 % of the bulge 
stars form in the late stage of evolution (a non-merger 
component) and also contribute to widen the distribution 
function. This is because the non-merger component of 
the bulge stars has much larger [Fe/H] than the merger 
component as explained below. 

The bulge stars can be divided into two chemically dif- 
ferent types with and without being contaminated by Type 
la supernova ejecta. With the adopted progenitor model 
(Kobayashi et al. 1998), the epoch of the first appearance 
of Type la supernovae tsNU = 0.5 Gyr after the forma- 
tion of each STAR. Thus, the STARs that formed after 
t ^ 0.5 Gyr are partially contaminated by Type la super- 
novae. Each Type la supernova ejects ~ 10 times more 
iron than a Type II supernova (Nomoto et al. 1997). As 
a result, iron contained in the non- merger component of 
the bulge stars partially originates from Type la super- 
novae and their [Fe/H] tends to be higher than that of the 
merger component. In Figure 12, the distribution function 
of the STARs that formed before and after t = 0.5 Gyr are 
plotted by the dashed and dotted lines, respectively, along 
with the total distribution function of the bulge stars (solid 
line). Both dashed and dotted lines show rather wide dis- 
tribution and clearly the peak value of each distribution 
is different. The very old STARs, which form before t = 
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0.5 Gyr, contribute mainly to the metal poor range of the 
distribution. 

The adopted yields are [0/Fe] = 0.4 from Type II super- 
novae and [0/Fe] = —1.5 from Type la supernovae (Table 
1). Thus the contamination by Type la supernovae ejecta 
decreases [0/Fe] in STARs starting from [0/Fe] = 0.4. In 
Figure 13, we present the metallicity distribution function 
of the bulge stars for [0/Fe]. The very peaky distribution 
shows that a half of the bulge stars have [0/Fe] ^ 0.4, i.e. 
a half of the bulge stars arc not contaminated by Type 
la supernovae. Also, almost 80 % of the bulge stars have 
[0/Fe] > 0. Figure 14 shows [0/Fe] as a function of [Fe/H] 
for the bulge stars. [Fe/H] increases up to ^ 0.3 without 
Type la supernova contamination. Later, due to the in- 
creasing Fe enrichment by Type la supernovae, [0/Fc] de- 
creases down to ^ —0.4 as [Fe/H] increases to ~ 0.7. We 
note for a given [Fc/H] (> —0.5), [0/Fe] is distributed over 
a wide range of > —0.4. Thus, even for [Fc/H] > 0, large 
fraction of the bulge stars have [0/Fo] > 0, while [0/Fc] 
is mostly sub-solar for [Fc/H] > 0.3. These features of 
our model are consistent with the observed results as plot- 
ted in Figure 14 by the dotted circles (Rich & McWilliam 
2000). The observed [Si/Fc] and [Ca/Fe] show also rela- 
tively large dispersion from —0.4 to -1-0.2 for [Fe/H] > 0, 
although the observed results are still preliminary (Rich & 
McWilliam 2000). 

In contrast, Matteucci & Brocato (1990) argued that 
the metallicity of the bulge stars mainly originated from 
Type II supernovae. In their model, they assumed tsma ~ 
1 Gyr and the star formation time scale is shorter than 
1 Gyr so that the majority of Type la supernovae have 
exploded long after the star formation. Therefore the one- 
zone chemical evolution model predicts that the Galactic 
bulge is a chemically single population. In contrast to this 
picture, our chemical and dynamical model predicts that 
the Galactic bulge is neither a chemically single population 
and a dynamically single population. 

5.3. Discussion 

There are several problems in our numerical model. 
First, the present model, which is compared with some 
observational data, is just the current best-fit model from 
only 5 realizations. Thus, the quantitative comparisons 
between our model and observation are still far from com- 
plete. One reason is that it is too difficult to find a consis- 
tent initial model for very good match to the observation 
before doing real calculations. Due to complicated merg- 
ing events occurring in the early phase, it is impossible 
to predict which initial halo will become a Milky- Way-like 
spiral galaxy. 

Another point is the resolution of the numerical model. 
As seen in Figure 7, the inner region of the present model 
is not well resolved. Thus, our definition radius for the 
bulge, R < 2 kpc, is not a clear definition when we com- 
pare our results with observational data. This is because 
the Baade's window, where the chemical properties of the 
bulge stars are observed, is located 0.5 kpc from the galac- 
tic center. The choice of our definition is determined from 
the limitation of the numerical resolution of the present 
model. We use the rather large radius of 2 kpc to define 
the bulge stars since the gravitational softening length of 
the present model is 0.5 kpc (in fact, the inspection of Fig- 



ure 7 shows that the effective resolution in our model is 
^ 1 kpc, which equals to the softening length of the dark 
matter particles). To summaries, the number of particles 
in the present model is still too small to fully resolve the 
chemical structure of the bulge. To make more sensible 
comparison between numerical models and the observa- 
tional data, we need much higher resolution models and 
more observational data obtained with large telescopes. 
Although there is such limitation, our model is in a good 
agreement with chemical and kinematic properties of the 
Galactic bulge stars. 

6. CONCLUSION 

In this paper, we report on the result of our high reso- 
lution three-dimensional models of the Galactic bulge for- 
mation. By modeling the formation and evolution of the 
Galaxy, we understand the formation history of the Galac- 
tic bulge as follows. 

The early evolution of the best-fit model reveals that 
most bulge stars form during the merger that occurred 
in the region of the deepest potential. Because of the 
strong star burst induced by the merger, the metallicity 
distribution of the bulge component becomes as wide as 
observed. Although such a strong star burst was a similar 
effect as the burst assumed in the previous chemical evo- 
lution model, our model show that there is another rea- 
son to widen the metallicity distribution function of the 
bulge. By analyzing our results, we find that the bulge 
stars in our model arc composed of two different compo- 
nents. One component (merger component) was produced 
by the merger occurred in the center of the halo as noted 
above. The other non-merger component was formed in 
the inner region of the disk after this merger. Due to the 
different formation epoch, the main source of iron for each 
component is different. Iron in the former comes mainly 
from Type II supernovae, and the other from Type la su- 
pernovae. Since each Type la supernova ejects ^ 10 times 
more iron than a Type II supernova, [Fe/H] of the non- 
merger component tends to be higher than that of the 
merger component, which widened the metallicity distri- 
bution function. Also, our model predicts a relatively large 
dispersion in [0/Fe] for [Fe/H] > 0. 

Although our model is not in a complete agreement with 
our Galaxy, some chemical and kinematic properties of the 
Galactic bulge, such as the velocity dispersion vs. metal- 
licity, are well reproduced. We thus conclude that an old 
fraction of the Galactic bulge is very likely to have been 
formed by the sub-clump merger in the proto-galactic en- 
vironment, while a younger fraction of the Galactic bulge 
has formed gradually in the inner disk. These two groups 
should show different chemical properties. 
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Table 1 

The iie.wv element ^melds lsed ik olh code 





M« 


Z 


Fe 





Type II SNe 


14.0 Mq 


2.54 Mq 


0.091 Mq 


1.80 Mq 


Type la SNe 


1.4 Mq 


1.4 Mq 


0.74 Mq 


0.14 Mq 



Note. — These data are taken from Tsujimoto et al. (1996) and Nomoto et al. (1997). 
"The total ejected mass per one supernova explosion. 



Table 2 

Kinematic properties of the model galaxies. 







spill (A) 


sicllar mass 


0.1 


0.59 


0.70 


4.2 X 10^" Mq 


0.5 


0.71 


0.52 


5.4 X 10^° Mq 


1.0 


0.77 


0.44 


4.6 X 10^° Mq 



Note. — ^ The star formation parameter. 
^ The ratio between the X and Z axis velocity dispersion. 
The stellar mass at f = 5 Gyr. 
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Fig. 1. The star formation probability (P) as a function of density of gas. The solid and dashed lines correspond to C = 1.0 and C 
0.1, respectively. 
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Fig. 2. — The star formation rate (SFR) as a function of time for galaxies with C = 0.1 (the dotted line) and C = 1.0 (the solid line). 
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Fig. 3. — The Type la supernova rate (ilsNia) as a function of time. 
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Fig. 4. — The schematic view of the feedback processes in our SPH code. "t(m)" is the stellar lifetime as a function of main-sequence mass 
(m/M©). 




Fig. 5. — The projected STAR positions at t = 5 Gyr (2 ~ 1) for the three components. From the left panel, the metal-poor STARs, the 
metal-rich and old STARs, and the young STARs are shown (see text). The size of the panels is 20kpc X 20kpc. 
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Fig. 6. — The star formation history for the bulge stars (solid line) and all stars (dotted line). The SFR divided by the total stellar mass 
of each population is shown as a function of time. 
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Fig. 7. — The face-on surface density of STARs at t = 5 Gyr. The crosses represent our model and the solid line is the fit to our model 
(see text). 




Fig. 8. — The projected STAR position for the first 0.5 Gyr of the evolution. During t = 0.4 - 0.5 Gyr, two clumps merge together. The 
size of the panels is lOkpc X lOkpc. 
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Fig. 9. — The mean velocity and the velocity dispersion for the bulge stars as a function of the distance from the galactic center. The upper 
and lower lines are the average velocity and the velocity dispersion, respectively. 
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Fig. 10. — The x-direction velocity dispersion for the bulge stars as a function [Fc/H]. The dotted circles are the observational data in 
Minniti (1996) (the error bound in the observation is ~ 8 - 14 km s~^). Our model STARs are shown by the filled squares (tform < 1 Gyr) 
and the filled triangles (tform < 5 Gyr). 
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Fig. 11. — The calculated mctallicity ([Fc/H]) distribution function of the bulge stars (solid line). The dashed line shows the observed 
metallicity distribution function of K-giant stars (McWilliam & Rich 1994). The distributions are normalized by its maximum value for 
compajrison. 
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Fig. 12. — The calculated metallicity ([Fe/H]) distribution function of the bulge stars (solid line). The distributions are not normalized. 
The dashed and dotted lines correspond to the very old STARs with tform < 0.5 Gyr and young STARs with tform > 0.5 Gyr, respectively. 
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Fig. 13. — The calculated metallicity ([O/Fe]) distribution function of the bulge stars. 
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Fig. 14. — [O/Fe] vs [Fe/H] for ~ 800 STARs, which are randomly selected from the bulge stars in our model (filled circles). The dotted 
circles show preliminary observational data (Rich & McWilliam 2000). 



